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We investigate bosonic Gaussian quantum states on an infinite cubic lattice in arbitrary spatial 



dimensions. We derive general properties of such states as ground states of quadratic Hamiltonians 
for both critical and non-critical cases. Tight analytic relations between the decay of the interaction 
O |' and the correlation functions are proven and the dependence of the correlation length on band 

gap and effective mass is derived. We show that properties of critical ground states depend on 
the gap of the point-symmetrized rather than on that of the original Hamiltonian. For critical 
systems with polynomially decaying interactions logarithmic deviations from polynomially decaying 
correlation functions are found. Moreover, we provide a generalization of the matrix product state 

■ representation for Gaussian states and show that properties hold analogously to the case of finite 

■ dimensional spin systems. 

\o ' 

I. INTRODUCTION 

0\ 

I The importance of bosonic Gaussian states arises from two facts. First, they provide a 

very good description for accessible states of a large variety of physical systems. In fact, 
every ground and thermal state of a quadratic bosonic Hamiltonian is Gaussian and remains 
so under quadratic time evolutions. In this way quadratic approximations naturally lead to 
Gaussian states. Hence, they are ubiquitous in quantum optics as well as in the description 
of vibrational modes in solid states, ion traps or nanomechanical oscillators. 

The second point for the relevance of Gaussian states is that they admit a powerful phase 
space description which enables us to solve quantum many-body problems which are other- 
wise (e.g., for spin systems) hardly tractable. In particular, the phase space dimension, and 
with it the complexity of many tasks, scales linearly rather than exponentially in the num- 
ber of involved subsystems. For this reason quadratic Hamiltonians and the corresponding 
' Gaussian states also play a paradigmatic role as they may serve as an exactly solvable toy 

model from which insight into other quantum systems may be gained. 

Exploiting the symplectic tools of the phase space description, exact solutions have been 
found for various problems in quantum information theory as well as in quantum statistical 
mechanics. In fact, many recent works form a bridge between these two fields as they address 
entanglement questions for asymptotically large lattices of quadratically coupled harmonic 
oscillators: the entropic area law 0. El has been investigated as well as entanglement 
statics [H 13- dynamics 0, H, @ and frustration 0, [n| . 

In the present paper we analytically derive general properties of ground states of trans- 
lationally invariant quadratic Hamiltonians on a cubic lattice. Moreover, we provide a 
representation of such states analogous to the matrix product states of finite dimensional 
spin systems. We start by giving an outlook and a non-technical summary of the main 
results. The results on the asymptotic scaling of ground state correlations are summarized 
in Tabled 

We note that related investigations of correlation functions were recently carried out in 
[T^ . fl3| for finite dimensional spin systems and in 0, ^| for generic harmonic lattices with 
non-critical finite range interactions. 
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TABLE I: Summary of the bounds derived in the paper on the asymptotic scaling of ground state 
correlations, depending on the scaling of the interaction (left column). Here n is the distance 
between two points (harmonic oscillators) on a cubic lattice of dimension d. O denotes upper 
bounds, O* tight upper bounds, and O the exact asyptotics. The table shows the results for generic 
interactions — special cases are discussed in the text. 

Quadratic Hamiltonians: In Sec. [HJ we start by introducing some basic results on 
quadratic Hamiltonians together with the used notation. 

TVanslationally invariant systems: In Sec. IIIII we show first that every pure transla- 
tional invariant Gaussian state is point symmetric. This implies that the spectral gap of the 
symmetrized rather than the original Hamiltonian determines the characteristic properties 
of the ground state. We provide a general formula for the latter and express its covariancc 
matrix in terms of a product of the inverse of the Fourier transformed spectral function and 
the Hamiltonian matrix. 

Non-critical systems: Sec. IIVI shows that if the Hamiltonian is gapped, then the correla- 
tions decay according to the interaction: a (super) polynomial decay of the interaction leads 
to the same (super) polynomial decay for the correlations, and (following Ref. Q) finite 
range interactions lead to exponentially decaying correlations. 

Correlation length and gap: Sec. fVl gives an explicit formula for the correlation length 
for gapped lD-Hamiltonians with finite range interactions. The correlation length £ is ex- 
pressed in terms of the dominating zero of the complex spectral function, which close to a 
critical point is in turn determined by the spectral gap A and the effective mass m* at the 
band gap via £ ~ (m*A) -1 / 2 . When the change in the Hamiltonian is given by a global 
scaling of the interactions this proves the folk theorem £ ~ 1/A. 

Critical systems: Sec. I VII shows that for generic <i-dimensional critical systems the corre- 
lations decay as 1/n +1 , where n is the distance between two points on the lattice. Whereas 
for sufficiently fast decreasing interactions in d = 1 the asymptotic bound is exactly poly- 
nomial, it contains an additional logarithmic correction for d > 2. Similarly for d = 1 a 
logarithmic deviation is found if the interaction decays exactly like —1/n 3 . 

Gaussian matrix product states: Sec. IVlll provides a representation of Gaussian states 
in terms of Gaussian matrix product states (GMPS). It is shown that any translational 
invariant pure state can be approximated by translational invariant GMPS to arbitrary pre- 
cision, that the correlations of any GMPS decay exponentially, and that every GMPS is the 
ground state of a local Hamiltonian. 

Appendix: Simulating Hamiltonians: We provide a Lemma showing that in d = 1 every 
translational invariant Hamiltonian can be simulated by any translational invariant nearest 
neighbor Hamiltonian supplemented by the set of translational invariant local Hamiltonians. 
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II. QUADRATIC HAMILTONIANS AND THEIR GROUND STATES 

Consider a system of N bosonic modes which are characterized by N pairs of canonical 
operators (Qi, Pi, ■ ■ ■ , Qn> Pn) —'■ R- The canonical commutation relations (CCR) are 
governed by the symplectic matrix a via 

[Rk, Rl] = i&kl , 

and the system may be equivalently described in terms of bosonic creation and annihilation 
operators a/ = (Qi + iP{)/^/2. Quadratic Hamiltonians are of the form 




T~L — - ^ HkiRkRi 



2 

kl 

where the Hamiltonian matrix H is real and positive semidefinite due to the Hermiticity and 
lower semi-boundedness of the Hamiltonian H.. Without loss of generality we neglect linear 
and constant terms since they can easily be incorporated by a displacement of the canonical 
operators and a change of the energy offset. Before we discuss the general case we mention 
some important special instances of quadratic Hamiltonians: a well studied ID example of 
this class is the case of nearest neighbor interactions in the position operators of harmonic 
oscillators on a chain with periodic boundary conditions 

1 N 

-H K = -Y,Q 2 i+P?-KQiQi+i> Ke[-i,i]. (i) 
»=i 

This kind of spring-like interaction was studied in the context of information transfer 0, 
entanglement statics and entanglement dynamics [jj. Moreover, it can be considered 

as the discretization of a massive bosonic continuum theory given by the Klein-Gordon 
Hamiltonian 



L/2 

' x) 2 + (V(j>(x)) Z + m 2 4>{x) 2 



L/2 



dx , 



where the coupling n is related to the mass m by = 1 + \{ II ^ ± ) 2 5}. Other finite 
range quadratic Hamiltonians appear as limiting cases of finite range spin Hamiltonians via 
the Holstein-Primakoff approximation [Tfij . In this way the xy-spin model with transverse 
magnetic field can for instance be mapped onto a quadratic bosonic Hamiltonian in the 
limit of strong polarization where a ~ {<r x + ia y )/2. Longer range interactions appear 
naturally for instance in ID systems of trapped ions. These can either be implemented as 
Coulomb crystals in Paul traps or in arrays of ion microtraps. When expanding around 
the equilibrium positions, the interaction between two ions at position i and j ^ i is — in 
harmonic approximation — of the form t^~h3 > where c > (c < 0) if Qi, Qj are position 
operators in radial (axial) direction 

Let us now return to the general case and briefly recall the normal mode decomposition 
[l7| : every Hamiltonian matrix can be brought to a diagonal normal form by a congruence 
transformation with a symplectic matrix S £ Sp(2iV, K) = {S\S(tS t = cr}: 1 



shs t = &) ^ u e(R , £l >o, (2) 




where the symplectic eigenvalues Si are the square roots of the duplicate nonzero eigenvalues 
of aHa T H . The diagonalizing symplectic transformation 5* has a unitary representation Us 



Note that we disregard systems where the Hamiltonian contains irrelevant normal modes. 
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on Hilbert space which transforms the Hamiltonian according to 

/ J i J 

ushuI = \ ^ {qi + p>) £l + i ]T p ] = E Ota + hh + 2 E ^ 2 • ( 3 ) 

Hence, by Eq. © the ground state energy E and the energy gap A can easily be expressed 
in terms of the symplectic eigenvalues of the Hamiltonian matrix: 

w - 1 c A - / min * Ei ' J = m 

^o-,^ £l , A-| 0j j>q . (4) 

The case of a vanishing energy gap A = is called critical and the respective ground states 
are often qualitatively different from those of non-critical Hamiltonians. For the Hamiltonian 
H K , Eq. Q), this happens in the strong coupling limit |k| = 1 — A 2 — > 1, and in the case 
of ID Coulomb crystals a vanishing energy gap in the radial modes can be considered as 
the origin of a structural phase transition where the linear alignment of the ions becomes 
unstable and changes to a zig-zag configuration 0,0; E3|- Needless to say that these phase 
transitions appear as well in higher dimensions and for various different configurations |2l| . 

Ground and thermal states of quadratic Hamiltonians arc Gaussian states, i.e, states hav- 
ing a Gaussian Wigner distribution in phas e sp ace. In the mathematical physics literature 
they are known as bosonic quasi-free states |22t l23j | . These states are completely character- 
ized by their first moments d k — tr [pi?fe] (which are w.l.o.g. set to zero in our case) and 
their covariance matrix (CM) 



Ikl = tr 



p{R k - d k ,Ri - di} , (5) 



where {•,•}+ is the anticommutator. The CM satisfies 7 > i<r, which expresses Heisenberg's 
uncertainty relation and is equivalent to the positivity of the corresponding density operator 
p > 0. In order to find the ground state of a quadratic Hamiltonian, observe that 



IE 



e t ® E = inf trlpH] @ \ inf tr[ 7 if] . (6) 

p 7 



By virtue of Eqs. (|2I3I) the infimum is attained for the ground state covariance matrix 

S, (7) 




lim S T 



which reduces to 7 = S T S in the non-critical case. Note that the ground state is unique as 
long as H does not contain irrelevant normal modes [which we have neglected from the very 
beginning in Eq. ©J ■ 

In many cases it is convenient to change the order of the canonical operators such that 
R = (Qii ■ ■ • j Qn, Pi, ■ ■ ■ j Pn)- Then the covariance matrix as well as the Hamiltonian 
matrix can be written in block form 

H = 

In this representation a quadratic Hamiltonian is particle number preserving iff Hq = Hp 
and Hqp = —Hq P , that is, the Hamiltonian contains only terms of the kind a\aj + a^ai. 

In quantum optics terms of the form a\a}-, which are not number preserving, are neglected 
within the framework of the rotating wave approximation. The resulting Hamiltonians have 
particular simple ground states: 




Theorem 1 (a) The ground state of any particle number preserving Hamiltonian is the 
vacuum with 7 = 1, and the corresponding ground state energy is given by Eq — jtrH . 
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Proof Number preserving Hamiltonians are most easily expressed in terms of creation and 
annihilation operators. For this reason we change to the respective complex representation 
via the transformation 




In this basis 77 is transformed to normal form via a block diagonal unitary transformation 
U ®U which in turn corresponds to an element of the orthogonal subgroup of the symplec- 
tic group Sp(27V,IR) n S0(2iV) ~ V(N) 24|. Hence, the diagonalizing S in Eqs. (00) is 
orthogonal and since J = due to particle number conservation, we have 7 = S T S = 1. Eq 
follows then immediately from Eq. 10. □ 

Another important class of quadratic Hamiltonians for which the ground state CM takes 
on a particular simple form corresponds to the case Hqp = where there is no coupling 
between the momentum and position operators: 

Theorem 1 (b) For a quadratic Hamiltonian with Hamiltonian matrix H — Hq © Hp the 
ground state energy and the ground state CM are given by 



1 = X®X-\ X = H Q 1 / 2 ^H 1 Q / 2 H P H 1 Q / 2 H- 1 ? 2 . (8) 



Proof Since aHa T H = HpHq © HqHp, the symplectic eigenvalues of H are given by the 
eigenvalues of w 'Hq\J 'Hp and thus Eq — ^tr [^/ Hg^/Hp] . Moreover, by the uniqueness of 
the ground state and the fact that Eq — jtr^i?] with 7 from Eq. (JSJ) we know that 7 is the 
ground state CM (as it is an admissible pure state CM by construction). □ 

Finally we give a general formula for the ground state CM in cases where the blocks in the 
Hamiltonian matrix can be diagonalized simultaneously. This is of particular importance as 
it applies to all translational invariant Hamiltonians discussed in the following sections. 

Theorem 1 (c) Consider a quadratic Hamiltonian for which the blocks Hq, Hp, Hqp of 
the Hamiltonian matrix can be diagonalized simultaneously and in addition Hqp = Hq P . 
Then with 

£ = ^HqHp — Hq P we have (9) 
E = itr[£] , A = X min (£) , 7 = {£ © £)- 1 oHo T . (10) 



Proof Since aHa T H — £ 2 © £ 2 we have indeed Eq = ^tv[£] and A = A m ; n (£'). Positivity 
7 > is implied by H > such that we can safely talk about the symplectic eigenvalues of 
7. The latter are, however, all equal to one due to (7c) 2 = —1 so that 7 is an admissible 
pure state CM. Moreover it belongs to the ground state since itr[i?7] = Eq. □ 



III. TRANSLATIONALLY INVARIANT SYSTEMS 



Let us now turn towards translationally invariant systems. We consider cubic lattices in 
d dimensions with periodic boundary conditions. For simplicity we assume that the size of 
the lattice is N d . The system is again characterized by a Hamiltonian matrix Hki, where 
the indices fc, I, which correspond to two points (harmonic oscillators) on the lattice, are 
now ci-component vectors in Zjy. Translational invariance is then reflected by the fact that 
any matrix element Am, A £ {Hq, Hp, Hqp} depends only on the relative position k — I 
of the two points on the lattice, and we will therefore often write Ak-i = A^i- Note that 
due to the periodic boundary conditions k — I is understood modulo N in each component. 
Matrices of this type are called circulant, and they are all simultaneously diagonalized via 
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the Fourier transform 



1 27Ti a 

T a g, = — = e N ap , a,BG T.n , such that 
VN 



A := jT® d AF t ® d = diag An e "^ m 

where mn is the usual scalar product in 7L d N . It follows immediately that all circulant 
matrices mutually commute. 

In the following, we will show that we can without loss of generality restrict ourselves 
to point-symmetric Hamiltonians, i.e., those for which Hqp — Hq P (which means that 
TC contains only pairs QuPi + QiPk)- Fot dimension d — 1 this is often called reflection 
symmetry. 

Theorem 2 Any translationally invariant pure state CM T is point symmetric. 



Proof For the proof, we use that any pure state covariance matrix can be written as 

r = ( r Q r QP \ = / x xy \ 
\rl P T P ) \yx x^ + yxy) ' 

where X > and Y is real and symmetric [25|. From translational invariance, it follows 
that all blocks and thus X and Y have to be circulant and therefore commute. Hence, 
Tqp = XY = YX = Tqp, i.e., T is point symmetric. □ 

Let V : Zfj Z d N be the reflection on the lattice and define the symmetrization operation 
S(A) = h(A + VAT) such that by the above theorem 5(7) = 7 for every translational 
invariant pure state CM. Then due to the cyclicity of the trace we have for any translational 
invariant Hamiltonian 

inftr[i7 7 l =inftr[5( J ff)7] . 

7 7 

Hence, the point-symmetrized Hamiltonian S(H), which differs from H by the off-diagonal 
block S(Hqp) = t;(Hqp + Hqp) has both the same ground state energy and the same 
ground state as H. Together with Theorem lc this leads us to the following: 

Theorem 3 Consider any translationally invariant quadratic Hamiltonian. With £ = 

1 /2 

[HqHp — j(Hqp + iJgp) 2 ] the ground state CM and the corresponding ground state 
energy are given by 

E = |tr[£] , 7 = (S 8 £)~ 1 aS(H)a T . (11) 

It is important to note that the energy gaps of H and S(H) will in general be different. In 
particular H might be gapless while S(H) is gapped. However, as we will see in the following 
sections, the properties of 7 depend on the gap A = A m i n (£) of the symmetrized Hamiltonian 
rather than on that of the original H. For this reason we will in the following for simplicity 
assume Hqp = Hq P . By Thm. 3 all results can then also be applied to the general case 
without point symmetry if one only keeps in mind that A is the gap corresponding to S(H). 

Note that the eigenvalues of £ are the symplectic eigenvalues oiS(H), i.e., £ — F® d £j?T® d 
is the excitation spectrum of the Hamiltonian. This is the reason for the notation where £ 
resides in Fourier space and £ in real space, which is differs from the normal usage of the 
hat throughout the paper. 

Correlation functions. According to Eqs. (|9I10I11I) we have to compute the entries of 
functions of matrices in order to learn about the entries of the covariance matrix i _This is most 
conveniently done by a double Fourier transformation, where one uses that f(M) = f(M), 
and we find 
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[f(M)] nm = -L^ e -^-[/(M)] rse ^ s ™ . (12) 



r. .s 



As we consider translationally invariant systems, M is circulant and thus M is diagonal. 
We define the function 

M(cj>) = M n e" m0 (13) 

such that M(2irr/N) = M r ^ r . As f(M) is solely determined by its first row, we can write 
[f(M)] n = ^- d J2 e 2mnr/N f(M(2nr/N)) . (14) 

In the following we will use the index n £ X d for the relative position of two points on 
the lattice. Their distance will be measured either by the or norm. Since we are 
considering finite dimensional lattices these are all equivalent for our purpose and we will 
simply write ||n||. In the thermodynamic limit N — > oo, the sum in Eq. (|14|l converges to 
the integral 

\J{M)] n = -l—[ d<!> f(M {(/>)) e in <t> with M{cj>) = V M n e~ in * , (15) 



(2ir) d 



T d — j 

1 r»GZ d 



where T d is the d-dimensional torus, i.e., [0,27r] d with periodic boundary conditions. The 
convergence holds as soon as \M n \ < oo [which holds e.g. for M n = 0(\\n\\~ a ) with some 
a > d] and / is continuous on an open interval which contains the range of M. 

From the definition (flTf of M, it follows that M g c S ?k (T d ) (the n times continuously 
differentiable functions on T d ) whenever the entries M n decay at least as fast as for 
some a > k + d, since then the sum of the derivatives converges uniformly. Particularly, if 
the entries of M decay faster than any polynomial, then M £ c S' oa {T d ). In the following the 
most important function of the type / o M will be the spectral function 




8 W = J L e_m0 IHqHpU - [E% P ]« • (16) 



Asymptotic notation. As the main issue of this paper is the asymptotic scaling of cor- 
relations, we use the Landau symbols o, O, and 0, as well as the symbol O* for tight 
bounds: 

• f(x) = o(g(x)) means lim ^j^l — 0, i.e., / vanishes strictly faster than g for x — * oo; 



f(x) — 0(g(x)), if lim sup 



fM 



is finite, i.e., / vanishes at least as fast as g; 



• f(x) = Q(g(x)), if f(x) = 0{g(x)) and g(x) — 0{f(x)) (i.e., exact asymptotics); 

• f(x) = 0*(g(x)), if f(x) = 0(g(xj) but f(x) ^ o(g(x)), i.e., g is a tight bound on 
f ? If / is taken from a set (e.g., those function consistent with the assumptions of a 
theorem) we will write / = 0*(g) if g is a tight bound for at least one / (i.e., the best 
possible universal bound under the given assumptions). 

If talking about Hamiltonians, the scaling is meant to hold for all blocks, e.g., if the in- 
teraction vanishes as 0(||n|| -Q ) for n — > oo, this holds for all the blocks Hq, Hp, and 
Hqp — HpQ. The same holds for covariance matrices in the non-critical case. By the 
shorthand notation f(n) — o(||n|| -00 ), we mean that f(n) = o(\\n\\~ a ) Va > 0. Note finally 
that the Landau symbols are also used in (Taylor) expansions around a point xq where the 
considered limit is x — > xq rather than x — > oo. 



2 In order to see the difference to ©, take an f(x) = g(x) for even x, f(x) = for odd x, x G N. Although 
/ does not bound g, thus f(x) ^ 0(g(x)), the bound g is certainly tight. A situation like this is met, e.g., 
in Theorem |5] where the correlations oscillate within an exponentially decaying envelope. 
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IV. NON-CRITICAL SYSTEMS 



In this section, we analyze the ground state correlations of non-critical systems, i.e., 
those which exhibit an energy gap A > between the ground an the first excited state. 
Simply speaking, we will show that the decay of correlations reflects the decay of the inter- 
action. While local (super-polynomially decaying) interactions imply exponentially (super- 
polynomially) decaying correlations, a polynomial decay of interactions will lead to the same 
polynomial law for the correlations. 

According to Theorem we will consider a translationally invariant system with a point- 
symmetric Hamiltonian (Hqp = Hq P ). Following I|1UI11|) . we have to determine the entries 

of (£- 1 8 £~ 1 )(jH<7 T ', with £ = (H Q H P + H^p) 1 / 2 . In Lemma H we will first show that 
it is possible to consider the two contributions independently, and as the asymptotics of 
ijH(t t is known, we only have to care about the entries of i.e., we have to determine 
the asymptotic behavior of the integral 

(^)« = ?Af / where £ = (H Q H P + H 2 QP ) 1/2 . 

Lemma 4 Given two asymptotic circulant matrices A, B in d dimensions with polynomially 
decaying entries, A n = 0(\\n\\~ a ), B n = 0(\\n\\~P), a, (3 > d. Then 

(AB) n = 0*(||n||-") , fi := mm{a,P} . 



Proof With Q v (n) := min{l, ||n|| n }, we know that \A n \ = 0(Q a ) and \B n \ — 0(Qp), and 



\(AB), 



3 



n-3. 



(17) 



We consider only one half space ||j|| < \\n — where we bound Qp(n — j) < Qp{n/2). As 
Qa(j) is summable, the contribution of this half-plane is 0{Qp{n/2)^. The other half-plane 
gives the same result with a and [3 interchanged, which proves the bound, while tightness 
follows by taking all A n , B n positive. □ 

We now determine the asymptotics of (£~ 1 ) n for different types of Hamiltonians. 

Lemma 5 For non-critical systems with rapidly decaying interactions, i.e., as o(|jn|| -00 ), 
the entries of ' £~ x decay rapidly as well. That is, 

A>0 (f- 1 ) n =o(||n||- 00 ) . 



Proof As the interactions decay as o(||7i|| -DO ), H. £ tf°°(T d ) (• = Q,P,PQ), and thus 
£ 2 = H Q H P + H 2 QP G (T d ). Since the system is gapped, i.e., £ > A > 0, it follows that 
also g := £~ 1 G c tf°°(T d ). For the proof, we need to bound 



(£-') 



(2tt) c 



d^g^e"^ 



by ||n|| K for all k £ N. First, let us have a look at the one-dimensional case. By integration 
by parts, we get 



Z7T 



i 



2nin 



d<t>g'(4>)e in * 



where the first part vanishes due to the periodicity of g. As g £ c ta ca (T 1 ), the integration 
by parts can be iterated arbitrarily often and all the brackets vanish, such that after k 
iterations, 



i 



27r(in) K 
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As g( K \<f>) is continuous, the integral can be bounded by / (</>)|d</> =: C K < oo, such 
that finally 

|(^ 1 )„|<% VkgN, 

which completes the proof of the one-dimensional case. 

The extension to higher dimensions is straightforward. For a given n = (nx, ■ ■ ■ ,nd), 
integrate by parts with respect to the fa for which \m\ = IMloo; we assume i = 1 without 
loss of generality. As g(-, <j>2, • ■ ■ , 4>d) € c ta°°(S 1 ) 1 the same arguments as in the ID case show 



(2tt)«*|th| 



QK 



c K 



□ 

For systems with local interactions, a stronger version of Lemma can be obtained: 
Lemma 6 For a system with finite range interaction, the entries of £~ l decay exponentially. 

This has been proven in Q for Hamiltonians of the type H = V © 1, exploiting a result 
on functions of banded matrices |26|. Following Eqs. I|9I11[1 the generalization to arbitrary 
translational invariant Hamiltonians is straightforward by replacing V with HqHp — Hq P . 
In fact, it has been shown recently that the result even extends to non translational invariant 
Hamiltonians of the form in Theorem 1 b [T3 |. 

Finally, we consider systems with polynomially decaying interaction. 

Lemma 7 For a ID lattice with H = V © 1 > and an exactly polynomially decaying 
interaction 

a b ,2<^N, 




E~ 1 decays polynomially with the same exponent, (£ _1 )„ = {V 1 ^ 2 ) n = ®{\n\~ v ). 

Hamiltonians of this type appear, e.g., for the vibrational degrees of freedom of ions in a 
linear trap, where v = 3. 

Proof We need to estimate (£ -1 )„ = (y~ 1/2 )„ = ^ V~ 1/2 (facade/) . Note that 

V r (^) = o + 26y"^^ = o + 26He[Li tf (e i *)l > , (18) 

71=1 

where Li„(z) = Xm>i z n jn v is the polylogarithm. The polynomial decay of coefficients 
implies V G ^"^(S 1 ), and as the system is non-critical, V^ 1/2 G C & V ~ 2 {S X ). As Li„ has 
an analytic continuation to C\[l; oo), V G & °°((0; 2tt)) and thus V' 1 ' 2 G ^°°((0; 2tt)). Wc 
can therefore integrate by parts v — 1 times, and as all brackets vanish due to periodicity, 
we obtain 



l 



2n(in) L 



2tt r 



—v- 1/2 (<t>) 



*64> , (19) 



and 



d^ 1 ^ 2y(0) 3 / 2 4^(0)5/2 

Note that the second term only appears if z/ > 3, and g only if > 4. As <?(0) G ^(S 1 ), 
its Fourier coefficients vanish as 0{n~ l ), as can be shown by integrating by parts. The 
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second term can be integrated by parts as well, the bracket vanishes due to continuity, and 
we remain with 



in 



with ft G "^(S 1 ). [For v = 3, a factor 2 appears as (F^)' = F^" 1 ).] As we will show later, 
y\y—i) j s absolutely integrable, hence the integral exists, and thus the Fourier coefficients 
of the second term in Eq. I|2U|) vanish as 0(n~ 1 ) as well. Finally, it remains to bound 



(21) 



/O 21/(0)3/ 2 

As Li^x) = Li„_i(:r)/iz:, it follows from Eq. I|18|) that 

V^Xcj)) = 26Re[r- 1 Li 1 (e 1 *)] = 2feRc [-i"' 1 log(l - e 4 *)] , 

where the last step is from the definition of Lii . 

We now distinguish two cases. First, assume that v is even. Then, 

oc Imlog(l - e 1 *) = arg(l - e^) = {<j>- tt)/2 

on (0;27r), hence the integrand in Eq. 121(1 is bounded and has a bounded derivative, and 
by integration by parts, the integral Eq. is 0(n _1 ). In case v is odd we have 

V*"- 1 ^) oc Relog(l-e # ) = log 1 1 - e**| = log(2 sin(0/2)) 
on (0; 27r). With h((f>) :— V~ 3 / 2 ((j))/2, the integrand in Eq. (|2"T)l can be written as 

V {v - 1] {(t>)h{4>) oc log(2 sin(</>/2)) h(0) + log(2 sin(</>/2)) [/i(0) - &(0)] . (22) 
The first term gives a contribution proportional to 

log(2 sin(0/2)) cos(n0)d(/) = — — 

2n 

as it is the back-transform of -|En>i cos W)/' 1 - ^or ^ ne secon d term, note that h £ 
^(S 1 ) for v > 3 and thus h(<f>) - h(Q) = ti {Q)<j> + o{4>) by Taylor's theorem. Therefore, the 
log singularity vanishes, and we can once more integrate by parts. The derivative is 

icot(0/2) - h(0)] +log(2sm(0/2)) ft'ty) . 

In the left part, the 1/0 singularity of cot(0/2) is cancelled out by h(4>) — h(Q) — 0(4>), and 
the second part is integrable as h! £ ^(S 1 ), so that the contribution of the integral (|21|l is 
0(n^ 1 ) as well. 

In order to show that n~ v is also a lower bound on (t/" 1 / 2 )^, one has to analyze the asymp- 
totics more carefully. Using the Riemann-Lebesgue lemma — which says that the Fourier co- 
efficients of absolutely integrable functions are o(l) — one finds that all terms in 119|) vanish 
as o(l/n 1/ ), except for the integral l|21l) . Now for even v, (|21|l can be integrated by parts, 
and while the brackets give a Q(n~ u ) term, the remaining integral is o(rt _l/ ), which proves 
that (V~ 1 / 2 ) n = 8(n~ l/ ). For odd v, on the other hand, the first part of (1221) gives exactly 
a polynomial decay, while the contributions from the second part vanishes as o(n~"), which 
proves {V^ 1 / 2 ) n = Q{rT lJ ) for odd v as well. □ 

Generalizations of Lemma [7J The preceding lemma can be extended to non-integer 
exponents a £" N: if V n oc n~ a , ti^O, then (£^ 1 ) n = 0(n~ a ). 

For the proof, define a = u + e, v £ N, < e < 1. Then V £ ^^(S 1 ), V £ ^°°((0; 2tt)), 
and one can integrate by parts v times, where all brackets vanish. What remains is to 
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bound the Fourier integral of the i/th derivative of V by n E . An upper bound can be 
established by noting that \V (l,) ((t>)\ < |Li e (e^)| = O^ 6-1 ) and \V^ +V >{(j))\ = 0(</> £ " 2 ). It 
follows that all contributions in the Fourier integral except the singularity from lead to 
o(l/n) contributions as can be shown by another integration by parts. In order to bound the 
Fourier integral of the 0(</> £_1 ) term, split the Fourier integral at ~. The integral over [0; i] 
can be directly bounded by n~ e , while for [^;1], an equivalent bound can be established 
after integration by parts, using 

y(u+X) = 0{<jf~ 2 ). This method is discussed in more detail 
in the proof of Theorem ^1 following Eq. I|44|) . 

The proof that n~ s is also a lower bound to (£~ 1 )„ is more involved. From a series 
expansion of V and its derivatives, it can be seen that it suffices to bound the sine and cosine 
Fourier coefficients of from below. As in the proof of Theorem 1 141 this is accomplished 
by splitting the integral into single oscillations of the sine or cosine and bounding each part 
by the derivative of <^ e_1 . 

For polynomially bounded interactions V n = 0(n~ a ), a > 1, not very much can be said 
without further knowledge. With v < a, v G N the largest integer strictly smaller than a, 
we know that V £ c (a l '~ 1 (S 1 ). Thus, one can integrate by parts v — 1 times, the brackets 
vanish, and the remaining Fourier integral is o(l) using the Riemann-Lebesgue lemma. It 
follows that (£~ 1 ) n — o(n~(" -1 )). In contrast to the case of an exactly polynomial decay, 
this can be extended to higher spatial dimensions d > 1 by replacing v—1 with v — d, which 
yields {E~ l ) n =o(n"^- d )). 

We now use the preceding lemmas about the entries of £ _1 (Lemma EH3 to derive corre- 
sponding results on the correlations of ground states of non-critical systems. 
Theorem 8 For systems with A > 0, the following holds: 

(i) If the Hamiltonian H has finite range, the ground state correlations decay exponen- 
tially. 

(ii) If H decays as o(||n|| _00 ) ; the ground state correlations decay as odl^H^ 00 ) as well. 

(Hi) For a ID system with H = V 1 where V decays with a power law v > 2, the 

ground state correlations decay as Q(\n\~ 1 '). 

Proof In all cases, we have to find the scaling of the ground state 7 which is the product 
7 = (£~ 1 (B£~ 1 )o-Ha T , Eq. (|1(J|) . Part (i) follows directly from Lemma|Hl as multiplying with 
a finite-range o~Ha T doesn't change the exponential decay, while (ii) follows from Lemma |S] 
the o(||n|| -00 ) decay of aHa T , and Lemma|| To show (Hi), note that for H = V ® 1, the 
ground state is 7 = 

y- 1 /2 yi/2 ; and from Lemma|3 0{n~ v ) follows. For f" 1 / 2 , LemmaQ 
also includes that the bound is exact, while for V 1 / 2 , it can be shown by transferring the 
proof of the lemma one-to-one. □ 

Note that a simple converse of Theorem [S] always holds: for each translationally invariant 
pure state CM 7, there exists a Hamiltonian H with the same asymptotic behavior as 7 
such that 7 is the ground state of H . This can be trivially seen by choosing H — a^/a T . 



V. CORRELATION LENGTH AND GAP 



In this section, we consider one-dimensional chains with local gapped Hamiltonians. We 
compute the correlation length for these systems and use this result to derive a relation 
between correlation length and gap. 



Theorem 9 Consider a non-critical ID chain with a local Hamiltonian. Define the complex 

\ L 1 1 / 2 

extension of the spectral function £((/)) — J2 n =a c ™ cos(n<j)) in Eq. flb\) as 



L 

£< 

n=0 
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FIG. 1: Sample arrangement of branch cuts and poles of ^fg inside the unit circle. From each odd 
order zero of g, a branch cut emerges. All cuts go to where they cancel with another cut. In case 
their number is odd, there is an additional branch point at cancelling the last cut. In case two 
zeros are on a line to the origin, the cuts are chosen curved. The integral of ^fg around the unit 
circle is equal to the integral around the cuts, plus integrals around the residues which originate 
from the even order zeros of g. 

such that g(e 1 ^) = £ 2 (<l>) = Hq{(f>)Hp{(j)) — Hq P ((/>) and let z be zero of g with the largest 
magnitude smaller than one. Then, the correlation length 

£ 1 



log |2| 

determines the asymptotic scaling of the correlations which is given by 

• O* (e~ n l^ I \pn), if z is a zero of order one, 

• 0*(e~ n ^), if z is a zero of even order, 

• o(e _ "/^ +e )) for all e > 0, if z is a zero of odd order larger than one. 

For the nearest neighbor interaction Hamiltonian H K from Eq. Q one has for instance 
£(4>) = yl — kcos(</>), so that g has simple zeros at zq = (1 ± vl — k 2 )/ k. Therefore 
5 = (1 — Vl — k 2 )/k, and the correlations decay as 0(e _ "^/i/n) where £ = — 1/ log \z\. 

Proof For local Hamiltonians, the correlations decay as the matrix elements of 
[Eq. Onj)]. By Fourier transforming @, S(<ft) = \fgJeW, with g(e^) = H Q {4>)H P {4>) - 
^Qpi't 1 ) = J2n=o c ™ cos ( n< t>) an even trigonometric polynomial (we assume cl 7^ without 
loss of generality), and min(<7(e 1 ^)) = A 2 . We have to compute 

1 f 27T 1 -a. 1 



2tt J £{(j)) 2m J s i y/g(z) 

where S 1 is the unit circle. The function g(z) has a pole of order L at zero and 2L zeros 
altogether. Since min(<?(0)) = A 2 > 0, g has no zeros on the unit circle. As g(z) = g(l/z), 
the zeros come in pairs, and L of them are inside the unit circle. Also, the conjugate of a 
zero is a zero as well. From each zero with odd multiplicity emerges a branch cut of \J g(z). 
We arrange all the branch cuts inside the unit circle such that they go straight to the middle 
where they annihilate with another cut. In case L is odd, the last cut is annihilated by 
the singularity of \/g(z) at 0. If two zeros lie on a line, one cut curves slightly. A sample 
arrangement is shown in Fig. ^ 

Following Cauchy's theorem, the integral can be decomposed into integrals along the 
different branch cuts and around the residues of 1/ ^fg, and one has to estimate the contri- 
butions from the different types of zeros of g. The simplest case is given by zeros zq with 
even multiplicity 2m. In that case, define h(z) := g(z)j{z— zq) 2 " 1 which has no zero around 
zq. The contribution from zq to the correlations is then given by the residue at zq and is 



d" 



(to - 1)! dz™- 1 \ y/hjz) 



(m—1) 

ocz 
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for n — (m — 1) > 0, i.e., it scales as \zo\ n . Note that for z$ $ R, the imaginary parts 
originating from zq and its conjugate z$ exactly cancel out, but the scaling is still given by 
\zq\ u = e nlog l 20 l, i.e., £ = — l/log|zo| is the corresponding correlation length. 

If zq is a simple zero of g(z), we have to integrate around the branch cut. Assume first 
that the cut goes to zero in a straight line, and consider a contour with distance e to the 
slit. Both the contribution from the e region around zero and the e semicircle at zq vanish 
as e — > 0, and the total integral is therefore given by twice the integral along the cut, 



1 



7™ JO ^Z - Z y/h(z) 



dz , 



where again h(z) = g(z)/(z — zq). Intuitively, for growing n the part of the integral close 
to zq becomes more and more dominating, i.e., the integral is well approximated by the 
modified integral where h(z) has been replaced by h(zo). After rotating it onto the real axis, 
this integral — up to a phase — reads 



1 



=dr 



- l ' 2 T{n) 



t^V\H z o)\ Jo ^/\z \ - r \A|/i(zo)| r(n + \) 



which for large n is 



y/n\zoh(za)\ V™ 



O 



IfoT 

„3/2 



(24) 



(25) 



In order to justify the approximation h(z) 
respective integrals. It is bounded by 



h(zo), consider the difference of the two 



1 



(*) 



On [zo/2,zo], h(z) is analytic and has no zeros, thus, |/i(z) -1 / 2 — h(zo)^ 1 / 2 \ < C\z — zq\, 
where C is the maximum of the derivative of h(z)" 1 ^ 2 on [zq/2, zq\. On [0, Zq/2], the same 
bound is obtained by choosing C the supremum of \h{z)^ x / 2 — h(zo)~ 1 / 2 \/\zo/2\ on [0, zq/2]. 
Together, (*) < C\z — zq\, and the above integral is bounded by 



C 



-M 



"-0 



rdr — C 



2T(n+|) 



O 



? 7 3/2 



i.e., it vanishes by l/n faster than the asymptotics derived in Eq. (|25|l . which justifies fixing 
h(z) at h(zo). 

From Eq. I|25|l . it follows that the scaling is e~ n ^ / ^/n, where the correlation length is 
again £ = — 1/ log | Zo \ ■ The same scaling behavior can be shown to hold for appropriately 
chosen curved branch cuts from z n to by relating the curved to a straight integral. 

The situation gets more complicated if zeros of odd order > 1 appear. In order to get an 
estimate which holds in all scenarios, we apply Cauchy's theorem to contract the unit circle 
in the integration (|23|) to a circle of radius r > \z \, where z is the largest zero inside the 
unit circle. Then, the integrand can be bounded by C r r n_1 (where C T < 00 is the supremum 
of lf-y/g on the circle), and this gives a bound 2-nC r r n ~ 1 for the integral. This holds for all 
r > I zo|, i.e., the correlation decay faster than e" losr for all r > |zq|. This does not imply 
that the correlations decay as e™ log ' Zo ', but it is still reasonable to define — 1 / log | zo | as the 
correlation length. □ 



Theorem 10 Consider a ID chain together with a family of Hamiltonians H(A) with gap 
A > 0, where H(A) is continuous for A — > in the sense that all entries of H converge. 
Then, the ground state correlations scale exponentially, and for sufficiently small A the 
correlation length is 
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Here, m* = ( 



is the effective mass at the band gap. 



For the discretized Klein-Gordon field JU, for example, we have A = ^/l — 1«|, 
m* = 2yjl — [k[/|k|, and for small A (corresponding to \k\ close to 1), one obtains 
£ ~ •v/]/t]/2(l — |/t]) ~ l/y/2A. Hence, the £ oc 1/A law holds if the coupling is increased 
relative to the on-site energy (in which case m* oc A). 

More generally, if we expand the spectral function [Eq. (|16|) ] around the band gap we are 
generically 3 led to the dispersion relation £(k) ~ yj A 2 + v 2 k 2 (k = <fi). By the definition of 
the effective mass and Theorem ^| this leads exactly to the folk theorem 

e-£. (26) 



Proof According to Theorem [§1 what remains to be done is to determine the position of 
the largest zero z of g in the unit circle. Due to the restriction on H(A), the coefficients 
of the polynomial g(z)z L and thus also the zeros of g continuously depend on A, i.e., for 
sufficiently small A, the zero closest to the unit circle is the one closest to the gap. In order 
to determine the position of this zero, we will expand g around the gap. We only discuss 
the generic case where the gap appears only for one angle </>o, g(4>o) = A. In the case of 
multiple occurrences of the gap in the spectrum, one will pick the gap which gives the zero 
closest to the unit circle, i.e., the largest correlation length. Furthermore, we assume <f>o = 
without loss of generality. Otherwise, one considers g(ze~ l ^°) instead of g(z), which on the 
unit circle coincides with the (rotated) spectrum. 

The knowledge on g =: u + iv (with u, v : C — ► R) which will be used in the proof is 

u(l) = A 2 «(1)=0 
u*(l) = v*(l)=0 (27) 

o^(l) = 2A/m* > «^(1)=0 

where the subscripts denote the partial derivative with respect to the respective subscript 
(in Euclidean coordinates z = x + iy, in polar coordinates z = re 1 ^). Note that z = 1 is the 
point where the gap appears, and that g(e 1 ^) = £((^) 2 is real. Therefore, the derivatives of 
the imaginary part v along the circle vanish, while the derivatives of the real part u are found 
to be u(l) = £(0) 2 = A 2 , u 4 (l) = 2£(0)£'(0) = 0, and u^(0) = 2£'(0) 2 + 2£(0)£"(0) = 
2 A/to*, where m* — 1 /£"((/>) is the effective mass at the band gap. 

We need to exploit the relation between Euclidean and polar coordinates, 

ffx(l)=Sr(l) ; 9y0-) =50(1) 

9xx(l) = g rr (l) ; #3,2,(1) =g^(l) +g r (l) 

and the Cauchy-Riemann equations u x — v y , u y = —v Xl and g xx + g yy = 0, which together 
with the information l|27l) lead to 

u(l) = A 2 ; v(l) = ; 
u x (l) = u y (l) = v x {\) = v y (l) = ; 
u xx {\) = -2A/m* ; u yy (l) = 2A/m* ; 
««:(!) = o ; %,(!)= o . 

Note that it is not possible to derive information about the mixed second derivates using 
only the information l|27|) . However, as long as v xy does not vanish at 1, v will only stay 



3 This makes the natural assumption that the minimum under the square root is quadratic. In fact, if it is of 
higher order, then m* = oo and thus £ = 0, which is consistent with the findings of the following section. 
An example of such a behavior is given by so called 'quadratic interactions' |jj for which H = V © 1, 
where V is the square of a banded matrix. 
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zero in direction of x or y, but not diagonally. Since A 2 > and 2A/m* > 0, the closest 
zero is — to second order — approximately located along the x axis. By intersecting with the 
parabola A 2 — -^{x — l) 2 , one finds that the zero is located at x$ w 1 — VAm*. For small 

A, the correlations thus decay with correlation length £ — l/log(l — VAmF) w l/yAm*. 

□ 



VI. CRITICAL SYSTEMS 

In the following, we discuss critical systems, i.e., systems without an energy gap, A = 
0. 4 In that case, the Hamiltonian will get singular and some entries of the ground state 
covariance matrix will diverge, which leads to difficulties and ambiguities in the description 
of the asymptotic behavior of correlations. We will therefore restrict to Hamiltonians of the 
type 

H=V®1, 

for which the ground state CM is 7 = V -1 / 2 © V 1 / 2 . While the Q part diverges, the entries 
of the P-block stay finite. Following Thm. 1(b) the extension to interactions of the form 
H = Hq © Hp is straight forward. 

In order to compute the correlations we have to determine the asymptotics of V 1 / 2 , i.e., 

We will restrict to the cases in which the excitation spectrum E = Vv" has only a finite 
number of zeros, i.e., finitely many points of criticality. In addition, we will also consider 
the special case in which the Hamiltonian exhibits a tensor product structure. 

We proceed as follows. First, we consider one-dimensional critical chains and show that 
the correlations decay typically as 0(n~ 2 ) and characterize those special cases where the 
correlations decay more rapidly. The practically important case of exactly cubic decaying 
interactions will be investigated in greater detail. Depending on the sign of the interaction 
this case will lead to a logarithmic deviation from the n~ 2 behavior. Then, we turn to higher 
dimensional systems and show that generically the correlations decay as logn, where 

d is the spatial dimension of the lattice. 

A. One dimension 

First, we prove a Lemma which shows that although taking the square root of a smooth 
function destroys its differentiability, the derivatives will stay bounded. 

Lemma 11 Let f g < ^" m ([— 1; 1]), f(x) > with the only zero at x — 0, and let 2u < m be 
the order of the minimum at x = 6, i.e., /( fc )(0) = Vfc < 2v, ,/ (2l/) (0) > 0. 
Define g(x) :— y / f(x). Then, the following holds: 

• For oddv, g g ^^([-l; 1]), and g g ^-"([-l; 0]), g g ^ m " I/ ([0; 1]), i.e., the first 
m — v derivatives (for x ^ ) are bounded. 

• For even v, g g ^"-"([-1; 1]). 



4 Note that there are different meanings of the notion criticality referring cither to a vanishing energy gap 
or to an algebraic decay of correlations. In this section we discuss in which cases these two properties are 
equivalent. 
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Proof Using the Taylor expansion f(x) = Y^k=2v °k xk + p( x )i ( x ) — o(x m k ) for k < to, 
we express g as g(x) = (sgn.x) v x l 'r(x) with 



r(x) 



\ 



k=2v 



CkX 



k-2v 



p{x) 



r 2u 



where we used that (sgn.x) l/ x l/ 



' . Let us now consider the derivatives of r(x). While 



the sum leads to a O(l) contribution, the fc'th derivative of the remainder behaves as 
o(l)/x 2u ~ m+k . Together, this leads to 



r ( fc )(x)=0(l) 2v 
r {k \x) = o{\) / x 2y -' m+k 2v 

Now consider the fc'th derivative of g(x) for x ^ 0, 



to + fc < , 
to + fc > 1 . 



g [k \x) = (sgnz) 



k 

E 

i=0 s. 



dx' 



Assume first k < v. Then, s ; oc 0(1)2^' for 2z/ — to + fc — Z < 0, and s ; cx o(l)a; m_ly_A; 
for 2^ — to + fc — Z > 1, and as to > 2v, it follows that </ fe ) — 0(a;) for fc < v, which 
cancels the discontinuity originating from sgnx. For k = u, on the contrary, Sk = 0(1), and 
sgnx introduces a discontinuity on g^ k \ yet, it remains bounded and piecewise differentiable 
on [— 1; 0] and [0; 1]. The first non-bounded si is found as soon as m — v — k = —1, and 
g g ^-"([O; 1]) directly follows. 

This also implies that for m— v— fc > 0, g(x)/(sgnx) 1 ' € c £" m ~ v 1]), i.e., the disconti- 
nuity is only due to (sgnx)^. Since, however, this is only discontinuous for odd v, it follows 
that g e <£ m - u {[-l] 1]) if v even. □ 



Theorem 12 Consider a one- dimensional critical chain with Hamiltonian H = V © 1, 
where V n = 0(n~ a ), a > 4 and where V has a finite number of critical points which are all 
quadratic minima ofV. Then, ( r )p) n = 0*{n~ 2 ). For V n oc n~ a , a > 3 it even follows that 
( 7 p) = 6(n- 2 ). 

Note that for V n oc n~ a , the extrema of V are always quadratic. 
Proof We want to estimate 

(V 1/2 ) n = y [ g^V^dt , (28) 

where g = V 1 ! 2 . Under both assumptions, V & < ^ 2 (5 1 ), and all critical points are minima of 
order 2. It follows from Lemma ITT1 that g is continuous with bounded derivative. Therefore, 
we can integrate by parts, the bracket vanishes, and we obtain 

(V 1/2 ) n = --±- r g\4>)e m ^4> . 

2mn Jo 

Now, split S 1 at the zeros of g into closed intervals Tj, [JjTj — S 1 , and rewrite the above 
integral as a sum of integrals over Tj, As g' <E c £{Tj) (and differentiable on the inner of Tj), 
one can once more integrate by parts which yields 

(Vl/2)n = ~2^fE ( - 9"(4>)e in *dtj . (29) 

Neither of the terms will vanish, but since g' £ c €(Tj\ the bracket is bounded. In case 
V n £ 0(n~ a ), a > 4, we have V £ ^ >3 ( l S 1 ), therefore g" is bounded (Lemma lllfl . and 
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the integrals vanish as o(l). Unless the contributions of the brackets for the different Ij 
cancel out, the n~~ 2 bound is tight, (V 1 ^ 2 ) n = 0*{rT 2 ). The tightness of the bound is also 
illustrated by the example which follows the proof. 

For the case of an exactly polynomial decay, we additionally have to show that g" is 
absolutely integrable for 3 < a < 4. Then, the exactness of the bound holds because the 
bracket in Eq. 1|29[) does not oscillate (the critical point is either at (j> = or at 4> = tt), 
and because the integral is o(l) for g" 6 In case the critical point is at (j> — tt, the 

latter holds since V G ^°°((0; 2tt)) implies that g" is bounded at tt, and V G ^(S 1 ) that 
g G < ^' 2 ((— 7r, tt)), which together proves that g" is bounded on S . 

In case the critical point is at (j) — 0, the situation is more involved (and for a = 3, 
a logarithmic correction appears, cf. Theorem I14|) . Since V^{4>) — — Im[Li Q _3(e 1 *)] = 
0((j> a ~ 4 ), we have 

V"{4>) = F"(0) + O(0 Q - 3 ), V\4>) = V"(0)<t> + O(0 Q - 2 ), V{(t>) = ±U"(O)0 2 + oi^- 1 ). 
With this information, 

which indeed proves that g" G C 1 ^ 1 ), and thus (V 1/2 ) n = 6(n~ 2 ). □ 

As an example, consider again the discretized Klein-Gordon field of Eq. which is critical 
for k = ±1, corresponding to V(4>) — 1 =F cos0. The Fourier integral is solvable and yields 
(7p)n = -^<^£ = Q(n- 2 ). 

Generalizations of Theorem 1121 Using Lemma ^TJ several generalizations for the ID 
critical case can be found. In the following, we mention some of them. In all cases H — V® 1 
is critical. 

Critical points of even order. — If V n = o(n~°°) and the critical points are minima of order 
2v, v even, the correlations decay as (jp) n = o(n~°°). This is the case, e.g., if V = X 2 with 
X itself rapidly decaying. 

Critical points of higher order. — If V has critical points of order at least 2v, v odd, and 
V n = 0(n- a ), a>2v + 2, then ( 7P )„ = 0( n -^ +1 '>). 

Minima of different orders. — If V has minima of different orders 2i/j, in general the mini- 
mum with the lowest odd ^ = vi will determine the asymptotics, (jp) n = Oin^^ 1 ^). As 
V S c € ly2 max l 1/ *}) (5 1 ) is required anyway, the piecewise differentiability of U 1 / 2 is guaran- 
teed. 

Weaker requirements on V. — It is possible to ease the requirements imposed on V in Thco- 
rem^Jto V n = 0(n~ a ), a > 3 or V G ^(S 1 ), respectively. The price one has to pay is that 
one gets an additional log correction as in the multidimensional critical case, Theorem 1151 
The method to bound g" is the same which is used there to derive (|39ll . 

The above proof does not cover the case of the relevant 1/n 3 interaction, which for instance 
appears for the motional degrees of freedom of trapped ions. In the following, we separately 
discuss this case. It will turn out that the scaling will depend on the sign of the coupling: 
while a positive sign (corresponding to the radial degrees of freedom) again gives a B(^j-) 
scaling as before, for the negative sign (corresponding to the axial degree of freedom) one 

getse(^P). 

Theorem 13 Consider a critical ID chain with a 1/n 3 coupling with positive sign, i.e., 
H = V ® t, V n — c/n 3 , Vo — 3cC(3)/2, c > 0, with £ the Riemann zeta function. Then, the 
ground state correlations scale as {^p) n = @(^2-)- 

Proof We take w.l.o.g. c = 1/2. For this sign of the coupling, the critical point is at 7r, 
V(tt) = 0. From the proof of Lemmad we know that V G ^{S 1 ), V G ^°°((0; 2tt)), and 
that V"{4>) = log(2sin(</>/2)) on (0;2tt). With g := U 1/2 , it follows from Lemma [TT] that 
g G ^(S 1 ), g G ^(K^])- and 9 € ^°°((0;tt]), g G <g°°{[— k; 0)). This means that all 
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derivatives g^ k \ k > 1 can exhibit jumps at the critical point n but they all remain bounded. 
In contrast, around 4> = 0, g' is continuous but g" has a log divergence. 
Thus, the Fourier integral 

(V 1/2 U = ±- I gWe^dcf, 
2tt J S i 

can be split at and 7r, and then integrated by parts twice. The brackets of the first 
integration cancel out due to continuity of g, and one remains with 

(V 1/2 ) n = (tftf) cos(n0)]J + J g"(<P) cos(n0)d</> 



where we used the symmetry of g. One finds [g'((j>) cos(n<^)]o = — y ^f^( — 1)™) an< ^ smce 
g" is integrable, the integral is o(l) due to the Riemann-Lebesgue lemma. Together, this 
proves (7p)„ = 0(4^). □ 

Theorem 14 Consider a critical ID chain with a l/n 3 coupling with negative sign, i.e., 
H = V © 1, V n = —c/n 3 , Vq = 2c£(3), c> 0, with £ the Riemann zeta function. Then, the 

ground state correlations scale as {~yp)n — 8( v/1 °f " -)■ 

Proof Again, take w.l.o.g. c = 1/2. For the negative sign of the interaction, the critical 
point is at = 0. Since at this point V" diverges, Lemma ITTI cannot be applied, and the 
situation gets more involved. 

As in the previous proof, we use that V S ^(S 1 ), V S ^°°((0; 2tt)), and thus V 1/2 S 
^(S 1 ), V 1 ' 2 € ^°°((0;27r)). Further, V"(<p) = - log(2 sin(0/2)) on (0;2tt), cf. the proof of 
Lcmma[71 and with sinx = x(l + 0(x 2 )) we have V"(4>) = — \og(<j>) + 0((j> 2 ) for (and 
similarly for <j> — > 27r) , and therefore 

V'(4>) = 0(l-log0)+O(0 3 ) , !/(0) = i0 2 (3-21og0) + O(0 4 ) . (30) 

As V 1 / 2 S ^(S 1 ), we can integrate by parts one time, 

(V 1/2 ) n = ^ / V- 1/2 {<pY n U<p = — r g'(cf>) sin(n0)d0 (31) 
2?r J 5 i 7rn J 

where we exploited the symmetry of V, and with g := V 1 / 2 . Then, from ()3()(l . 



and after another round of approximation, 



= + Of ) , g'W = -^-^= + O 



>/2 V^loi^y ' ' 23/2 ^yfl^l lo g 0| 3 / 2 



This shows that the remainder <?'(</>) — yj log </>|/2 is continuous with a absolutely integrable 
derivative, and by integration by parts it follows that it only leads to a contribution 0(l/n) 
in the integral l|31|) . Thus, it remains to investigate the asymptotics of the sine Fourier 
coefficients of h((f>) — y/\ log^>|. For convenience, we split the integral l|31|) at 1, and [1;tt] 
only contributes with 0(l/n), as h is continuous with absolutely integrable derivative on 
[l;7r]. On [0; 1], we have to compute the asymptotics of 



1=1 \/ -\og<l>sm{n(j))A4> . (32) 
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Therefore, split the integral at 1/n. The left integral can be bounded directly, and the right 
after integration by parts [cf. the treatment of Eq. (|44H ]. One gets 



X < 



V-log^d0+ + - / — d0 = O 



In order to prove that this is also a lower bound for the asymptotics, it suffices to show this 
for the integral as all other contributions vanish more quickly. To this end, split the 
integral (|32[1 into single oscillations of the sine, Jk = [^2*, 2?r (fc+ 1 ) > o. As -y/— log has 
negative slope on (0; 1), each of the Jk gives a positive contribution to X, and thus we can 
truncate the integral at i, 



X > ^2 / V- lo S0 sin 



(33) 



On [0; 5], \/ — log has a positive curvature, and thus, each of the integrals can be estimated 
by linearly approximating y/— log cj) at the middle of each Jk but with the slope at 27r (k+i) , 
which gives 

V ' — log0 sin(n0) d</> > — j 



2 

Jk 71 2^(fc+l) 



-log 



2^r(fc+l) 



2 



Now, we plug this into the sum (|33|) and bound the sum by the integral from 2zl to 
(the integrand in monotonically decreasing), which indeed gives a lower bound ^{yjlog 3^ 
\/log 2) on X and thus proves the 9(\/logn/n 2 ) scaling. □ 



B. Higher dimensions 

For more than one dimension, the situation is more involved. First of all, it is clear by 
taking many uncoupled copies of the one-dimensional chain that there exist cases where the 
correlations will decay as n~ 2 . However, these are very special examples corresponding to 
Hamiltonians with a tensor product structure Hi 1 i 2 ^ 1 j 2 — Hi 1 j 1 H' i2 - 2 . In contrast, we show 
that for generic systems the correlations in the critical case decay as 0{n~^ d+1 ^ logn), where 
d is the dimension of the lattice. The requirement is again that the energy spectrum £(</>) 
has only a finite number of zeroes, i.e., finitely many critical points. 

Note that the case of a Hamiltonian with a tensor product structure can also be solved, 
as in that case V becomes a product of terms depending on one <pi each and thus the 
integral factorizes. Interestingly, although the correlations along the axes decay as n~ 2 , the 
correlations in a fixed diagonal direction will decay as n^ 2 ■ ■ ■ n^ 2 oc ||n||~ 2<i and thus even 
faster than in the following theorem. The 0(|jn||~( d+1 ) log ||n||) decay of the theorem holds 
isotropically, i.e., independent of the direction of n. 

Theorem 15 Consider a d- dimensional bosonic lattice with a critical Hamiltonian H = 
V © A . Then the P- correlations of the ground state decay as 



O 



(j|n||-( d+1 >log||n|| 



if the following holds: V G < r£ d+1 [e.g., the correlations decay as 0(||n|| ( 2d + 1 + £ )) ; e>Q],V 
has only a finite number of zeros which are quadratic minima, i.e., the Hessian y g^g^ ^ 
is positive definite at all zeros. 

Proof We have to evaluate the asymptotic behavior of the integral 
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Let us first briefly sketch the proof. We start by showing that it suffices to analyze each 
critical point separately. To this end, we show that is is possible to smoothly cut out some 
environment of each critical point which reproduces the asymptotic behavior. Then, we 
rotate the coordinate system such that we always look at the correlations in a fixed direction, 
and integrate by parts — which surprisingly can be carried out as often as V is differentiable, 
as all the brackets vanish. Therefore, the information about the asymptotics is contained in 
the remaining integral, and after a properly chosen number of partial integrations, we will 
attempt to estimate this term. 

Let now d, i = 1, . . . , J be the zeros of V. Clearly, these will be the only points which 

contribute to the asymptotics as everywhere else \fv is c tc? d+1 . In order to separate the 
contributions coming from the different £,-, we will make use of so-called neutralizers |27| . 
For our purposes, these are functions Af^ 0>r S ^°°(M. d — > [0; 1]) which satifsy 



1 : 11? -Coll <r/2 
0:||£-£o||>r 



and are rotationally symmetric (cf. 27] for an explicit construction). For each Q, there 
exists an r, such that the balls B ri (Q) do not intersect. We now define the functions 

Clearly, p is c ta d+1 , and so is each fa except at Q. Furthermore, each fi is still the square 
root of a 'rf d+1 function. By definition, 

0^ 1/2 )» = 7AdE / ^(^cosM + tAi / d d #(0)cosM, (34) 

(Z7T) . =l J T d (Z7TJ J T d 

i.e., it suffices to look at the asymptotics of each fi separately. The contribution of p is 
0(||n|| - ( d+1 )) as can be shown by successive integrations by parts just as for the non-critial 
lattice (cf. the proof of Lemma [SJ). 
Let us now analyze the integrals 

d^^cosM. 

B n (Ci) 

The integration range can be restriced to B ri (Q) as /j vanishes outside the ball. By a 
rotation, this can be mapped to an integral where n = (\\n\\, 0, . . . , 0), whereas fi is rotated 
to another function /, with the same properties, 

h= I d d 0A(0) cos[||»||0i] . 

Since the integrand is continuous and thus bounded, it is absolutely integrable, and from 
Fubini's theorem, one finds 

li = J d d ^4> 

BrAGi) C '' 1 " r - 



where we separated out the integration over the first component. The vector <f> denotes the 
components 2 ... d of <j>. The extension of the integration range to a cylinder is possible as 
/, vanishes outside B n (£,-). 

Let us now require (j> ^ Q. This does not change the integral since the excluded set is of 
measure zero, but it ensures that fi is in < W d+1 . This allows us to integrate the inner integral 
Ji{4>) by parts up to d + 1 times, and each of the brackets 



/^(^^^(IM!^ - W2) 



Ci.i+r-i 

01=Ci,l-»'i 
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appearing in the fc'th integration step vanishes. Here, f} (<f>i,4>) = d d fi(4>\, <f>)/d<pf is the 
d'th partial derivative with respect to the first argument. After integrating by parts d times, 
we obtain 



0,i + n 

U f AM. tW 



B Ti (Ci) Ci,i 



fl a > (fa, 0)cos[||n||^-d7r/2] . (35) 



Now we proceed as follows: first, we show that the order of integration can be interchanged, 
and second, we show that for the function obtained after integrating over 0, the Fourier 
coefficients vanish as log(||n||)/||n||. 

The central issue for what follows is to find suitable bounds on \fi \- Therefore, define 
ff =: hi £ c ta d+1 . By virtue of Taylor's theorem, and as hi(d) = is a minimum, 

hi(<t>) = U ( t>- d) ■ (D a k(Ci))(0 - Ci) + o(H0 - on 2 ) 

with D 2 the second derivative. As the first term is bounded by ^ ||D 2 /ii (Ci) II oo II — Oil 2 an d 
the second vanished faster than \\cj> — d\\ 2 , we can find E; L > and C\ > such that 

Mfal^dU-Cif V\\<f>-Ci\\<ei. (36) 

By looking at the Tayor series of h\ = dhi/dfa up to the first order we also find that there 
are > and C 2 > such that 

\h'M\ <C 2 \\<f>- Ci|| V||0 — C*ll < e* - (37) 

In addition to these upper bounds, we will also need a lower bound on \hi\. Again, by the 
Taylor expansion of hi around Q , we find 

\hi{ct>)\ > A min [B 2 h t (Q)]~o(U~Q\\ 2 ) , 

and as all the zeros are quadratic minima, i.e., A m i n [D 2 /ii(Ci)] > 0, there exist Si > 0, 
C3 > such that 

M4>)\ > C 3 ||0 — C.ll 2 V||0-Ci|| <ei . (38) 

Clearly, e$ can be chosen equal in Eqs. I|3t)l38|l . Note that the bounds can be chosen to 
be invariant under rotation of hi and thus of /j. This holds in particular for the as 
the remainders of Taylor series vanish uniformly. Thus, the bound we will obtain for the 
correlation function indeed only depends on ||n|| and not on the direction of n. 

Now, we use the conditions (|36I38|) to derive bounds on \ \. Therefore, note that from 
fi = \fhl it follows that 

E c n .., k h^---h^ 

JlH Hk= k 

7(k) _ ji/=o,i, 2,... 
J i 



,(2fe-l)/2 

One can easily check that for each term in the numerator, the number Kq of zeroth deriva- 
tives and the number K\ of first derivatives of hi satisfy 2Kq + K\ > k. By bounding all 
higher derivatives of hi from above by constants, we find that the modulus of each summand 
in the numerator, and thus the modulus of the numerator itself, can be bounded above by 
C"||</>— Ci\\ k in the ball B ei (Q) with some C > 0. On the other hand, it follows directly from 
P5|l that the modulus of the denominator is bounded below by C"\\<j> — Cj|| 2fe- \ C" > 0, 
such that in total 

\f?\<l>)\<C 1 k-1 ; l<A<d + l. (39) 



Note that this holds not only inside B €i (d) but in the whole domain of fi, as outside B ei (Q), 
fi is ( & d+1 and thus all the derivatives are bounded. 
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Eq. I|39[) is the key result for the remaining part of the proof. First, it can be used to 
bound the integrand in (|35|l by an integrable singularity (this is most easily seen in spherical 
coordinates, where l/r d_1 is integrable in a d-dimensional space). Hence, the order of 
integration in H35[) can be interchanged, and it remains to investigate the asymptotics of the 
integral 



1 



0,1 + r i 



HI". 

0,i - n 



i#i(</>i)cos[||n|l<?(>i - dn/2] , with 



(40) 
(41) 



From l|39(l . we now derive bounds on gi(4>i) and its first derivative. Again, we may safely 
fix ^ Q i as this has measure zero. Then, using l|39[) we find that 



WVPi, 



< 



C 



lo ((0i-Cm) 2 +^ 2 )) ( ^ 1)/2 

where we have transformed into spherical coordinates [Sd-i is the surface of the (d — 1)- 
dimensional unit sphere] and assumed the Z 2 -norm. Since {4>\ — Ci) 2 +r 2 > r 2 , the integrand 
can be bounded once again, and we find 

CSd-i 



IS* W>i ) I < 



o ((0i-Co) 2 +r 2 ) 1/2 
- C\-\og\cf> l -C i ,i\+ log 
< -Clogl0i-G.il 



dr 



Cm) 5 



(42) 



where in the last step we used that in (|40|) |0i — Cill < r % an d that can be chosen 
sufficiently small. 

Next, we derive a bound on g[((f)i). As we fix <f>\ ^ Ci) the integrand in l|41|l is ^ and we 
can take the differentiation into the integral, 



g'M 



Again, we bound the integrand by virtue of Eq. (|39() and obtain 

CSd-i 



IftOMI < 



o ((0i-O,i) 2 + ^ 2 ) 
arctan 



dr 



C- 



I0l-Ci,l| 



< 



c 



i^i - Ci.ii 101 - c»,ii ' 43 

Finally, these two bounds will allow us to estimate (|40|) and thus the asymptotics of the 
correlations in the lattice. We consider one half of the integral H40p . 

0,1 + r i 

d0i5i(0i)cos[||n||0i-d7r/2] , (44) 



0,i 

as both halves contribute equally to the asymptotics. We then split the integral at i 
r,-/||n||. The left part gives 



Ci,i+n/||n|| 



d0igi(0i)cos[||n||0i -dn/2] 



EH 

< 



l(— log |01 - 0,1 1) 



Ci,i 



Tj - Tj \ogr l + Tj log [jnjj 

Hi 



(45) 
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The right part of the split integral (|44|l can be estimated by integration by parts, 



d(/)i3i(^i)cos[||n||(/)i -d-Tt/2} 



< 



< 



5i(^i)-n — rr cos[||tt.||0i - (d + 1)71-/2] 



1 d.i+n 



Ci,l+Ti/||n| 



\\n\\ 

Ci,l+ r i/IMI 



d4>M<h)\ 



™3 c log| 



(46) 



Thus, both halves [Eqs. (051,103] S ive a l°g IMI/ll n ll bound for the integral 03J), and thus 
the integral Ii is asymptotically bounded by log ||«-|l/IM| d+1 following Eq. (|4Ufl. As the 
number of such integrals in (fTT^f) is finite, this proves that the correlations of the ground 
state decay at least as log ||n||/||n|| d+1 . □ 



VII. GAUSSIAN MATRIX PRODUCT STATES 



Recently, so-called Matrix Product States (MPS) have attracted growing interest in quan- 
tum information theory. These states appear in the DMRG (Density Matrix Renormaliza- 
tion Group) method which is a powerful tool to compute ground state properties of transla- 
tional invariant Hamiltonians. From a quantum information perspective, this class of states 
can be given a physical interpretation in terms of projected entangled pairs: they can be 
obtained by taking a chain of maximally entangled pairs of dimension D and applying a 
map T as in Fig. [21 in a translational invariant fashion. In the limit of large bond dimension 
D, this allows to approximate arbitrary translational invariant states. In finite dimensions, 
the MPS representation turned out to be a very fruitful concept as it led to new power- 
ful numerical algorithms [2^, 12^. I30L l3lj accompanied by a better understan ding of their 
efficiency |32| . and new insights into renormalization group transformations |33j and se- 
quential quantum generators |34|. In the following, we generalize matrix product states to 
the Gaussian scenario. 



A. Definition of Gaussian MPS 



We start by defining Gaussian matrix product states (GMPS). The definition resembles 
the physical interpretation of finite-dimensional matrix product states as projected entangled 
pairs (PEPs). In finite dimensions, MPS can be described by taking maximally entangled 
pairs of dimension D between adjacent sites, and applying arbitrary local operations on each 
site, mapping the DxD input to a <i-dimensional output state. Similarly, GMPS are obtained 
by taking a number of entangled bonds and applying local (not necessarily trace-preserving) 




FIG. 2: Construction of Gaussian Matrix Product States (GMPS). GMPS are obtained by taking 
a fixed number M of maximally entangled (i.e., EPR) pairs shared by adjacent sites, and applying 
an arbitrary 2M to 1 mode Gaussian operation T' ! ' on site i. 
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Yout 

FIG. 3: The Jamiolkowski isomorphism. The Gaussian channel described by the state F can be 
implemented by projecting the input state 7i n (mode A) and the input port of F (mode B) onto 
the EPR state (symbolized by curly brackets). In case of success, the output is obtained in mode 
C. The operation can be made trace-preserving by measuring in a basis of displaced EPR states, 
and displacing C according to the measurement outcome. 



operations , where the boundary conditions can be taken either open or closed. Any 
GMPS is completely described by the type of the bonds and by the operations 7~M. Note 
that this construction holds independent of the spatial dimension. For one dimension, it is 
illustrated in Fig. [21 As matrix product states are frequently used to describe translationally 
invariant systems, an inportant case is given if all maps are identical, 7"W = T Vi. 

In order to define MPS in the Gaussian world, we have to decide on the type of the bonds 
as well as on the type of operations. We choose both the bonds to be Gaussian states and the 
operations to be Gaussian operations, i.e., operations mapping Gaussian inputs to Gaussian 
outputs. For now, we will take the bonds to be maximally entangled (i.e., EPR) states, such 
that the only parameter originating from the bonds is the number M of EPRs. We show 
later on how the case of finitely entangled bonds can be easily embedded. 

As to the operations, we will allow for arbitrary Gaussian operations. Operations of this 
type are most easily described by the Jamiolkowski isomorphism [35|. There, any Gaussian 
operation T which maps N input modes to M output modes can be described by an N + M 
mode covariance matrix T with block B (input) and C (output). The corresponding map 
on some input state 7; n in mode A is implemented by projecting the modes A and B onto 
an EPR state as shown in Fig. El such that the output state T(y ln ) is obtained in mode C. 
Conversely, the matrix Y which represents the channel T is obtained by applying the channel 
to one half of a maximally entangled state. The duality between T and T is most easily 
understood in terms of teleportation, and shows that this characterization encompasses all 
Gaussian operations. Note that the protocol of Fig.|3]can be always made trace-preserving by 
projecting onto the set of phase-space displaced EPR states and correcting the displacement 
of mode C according to the measurement outcome |86j. 

In the following, we will denote all maps T by their corresponding CM T. Sometimes, we 
will speak of the modes B and C as input and output ports of T, respectively. 

We now discuss how the covariance matrix of the output will depend on the CM of the 
input and on the channel T |36l l37j . This is most easily computed in the framework of 
characteristic functions [3^| . The characteristic function of the output is given by 

Xc(£c)« J e-^ u e- e *° v te c 8{x A -x B )8{pA+PB)&UB , 
and by integrating out subsystem A, 

Xc(£c)<x / 'e-^ C MfB Cd £ B , 

with 

M= ( Ojo + r B t bc \ 

Basically, the integration J c\^aS(xa — xb)5{pa + Pb) does the following: first, it applies 
the partial transposition 6 = ( o -i ) to one of the subsystems, and second, it collapses the 
two systems A and B in the covariance matrix by adding the corresponding entries. The 
integration over £b, one the other hand, leads to a state whose CM is the Schur complement 
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of Mn, M 2 2 - M 2 iM 1 " 1 1 Mi2, such that the output state is described by the CM 

Tout = r c — T CB - — — — „F BC . 

Let us briefly summarize how to perform projective measurements onto the EPR state in 
the framework of CMs, where we denote the measured modes by A and B, while C is the 
remaining part of the system. First, apply the partial transposition to B, second, collapse 
A and B, and third, take the Schur complement of the collapsed mode AB, which gives the 
output CM of C. 

As we discuss Gaussian matrix product states in connection with ground states of Hamil- 
tonians, we are mainly interested in pure GMPS. Particularly, a GMPS is pure if the I™ 
which describe the operations are taken to be pure, which we assume from now on. 

Let us finally emphasize that the given defintion of MPS holds independent of the spatial 
dimension of the system, as do most of the following results, and in fact applies to an 
arbitrary graph. 



B. Completeness of Gaussian MPS 



In the following, we show that any pure and translational invariant state can be approxi- 
mated arbitrarily well by translational invariant Gaussian matrix product states, i.e., GMPS 
with identical local operations T. (Without translational invariance, this is clear anyway: 
the complete state is prepared locally and teleported to its destination using the bonds.) The 
proof is presented for one dimension, but can be extended to higher spatial dimensions. For 
the proof, we use a theorem on the simulation of translational invariant Hamiltonians which 
is proven in the Appendix (Thcorcm llol . For our purposes, it says that in order to simulate 
an arbitrary translational invariant Hamiltionian with reflection symmetry, it suffices if one 
can implement translational invariant local and nearest neighbor Hamiltonians. 5 




FIG. 4: Implementation of a translational invariant nearest neighbor Hamiltonian in a translational 
invariant fashion. Starting from 7i n , the input ist first teleported to the left, then, the infinitesimal 
time evolution S = e aH , H <C 1, is performed, and finally, the state is teleported back. 



5 This naturally extends to higher dimensions. For two dimensions, e.g., one needs nearest neighbor inter- 
actions and in addition one interaction with the closest neighbor along the diagonal in order to break the 
reflection symmetry. The implementation of such an interaction in the MPS formalism is a straightforward 
extension of the presented method. 



26 



Given a translational invariant state 7, there is a translational invariant Hamiltonian 
H which transforms the separable state 1 into 7, 7 = SS T , S — e aH . According to 
Theorem 1161 this time evolution can be approximated arbitrarily well by a sequence of 
translational invariant local (one-mode) and nearest neighbor (two-mode) Hamiltonians Hj , 

e'^jje®.^, (47) 

where the Hj act on one or two modes, respectively, and approach the identity for growing 
J. 

Clearly, translational invariant local Hamiltionians can be implemented by local maps 
without using any EPR bonds. In the following, we show how translational invariant nearest- 
neighbor interactions can be implemented by exploiting the entanglement of the bonds. The 
whole procedure is illustrated in Fig. 0] and requires two EPR pairs per site. We start with 
some initial state 7i n onto which we want to apply 5® = e® aHj w 1 + ® aHj . 

First, we perform local EPR measurements between the modes of 7i n and one of the bonds 
in order to teleport the modes of 7; n to the left, cf. Fig.^J,. Then, the infinitesimal symplectic 
operation S = e aHj is applied to the left-teleported mode and the second bond, Fig. 0Jd. 
In the last step, another EPR measurement is performed which teleports the left-tclcportcd 
mode back to the right, and "into" the mode on which the adjacent S was applied. As the 
operations e aHj « 1 + aHj all commute, the "nested" application of the nearest neighbor 
symplectic operations S indeed give iS®. Thus, the remaining mode indeed contains the 
output 7 out = S'®7i n S , 0. The whole decomposition (|47|l can be implemented by iterated 
application of the whole protocol of Fig. 

C. Gaussian MPS with finitely entangled bonds 

In this subsection we show that in general, infinitely entangled bonds can be replaced by 
finitely entangled ones. Intuitively, this should be possible whenever the channel destroys 
some of the entanglement of the bond anyway, i.e., rW is non-maximally entangled. In that 
case, it should be possible to use a less entangled bond while choosing a channel which does 
not destroy entanglement any more. 

The method is illustrated in Fig. |SJ Again, for reasons of clarity we restrict to one 
dimension and one bond. The argument however appies independent of the spatial dimension 



i+l 




FIG. 5: How to make the bonds of GMPS finitely entangled, a) The initial MPS. b) Do a Schmidt 
decomposition of the original map F. c) Move the Sj^ through the infinitely entangled bond to the 
next site, d) Swap the finitely and the infinitely entangled pair. 
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and the number of bonds. The only restriction we have to make is the restriction to pure 
GMPS, i.e., those with pure T®. 

Consider a GMPS with local channe ls g iven by and infinitely entangled bonds, Fig.^. 
First, apply a Schmidt decomposition 39] to rH in the partition A\BC, which can be always 
done as long as 1™ is pure. The Schmidt decomposition allows us to rewrite the state as 
shown in Fig. [SJj — an entangled state between modes A and C with two-mode squeezing 
rW, B in the coherent state 1, and sympectic operations and S^ c which are applied to 
modes A and BC, respectively. As the bond itself is infinitely entangled, we can teleport 
the sympectic operation through the bond to the next site as indicated in Fig. [SJd. Then, 
can b e mer g e d with S^ c to a new operation acting on modes B and C of site 
i (Fig. IHt)- Finally, in the triples consisting of one maximally entangled state, one non- 
maximally entangled state, and the projection onto the EPR state, the maximally and the 
non-maximally entangled state can be swapped, resulting in Fig. [SJi. There, we have finitely 
entangled bonds, while the infinite entanglement has been moved into the new maps rM. 

It is tempting to apply this construction to the completeness proof of the preceding section 
in order to obtain a construction which is less wasting with respect to resources. However, 
for any iterative protocol this is most likely difficult to achieve. The reason for this is 
found in the no-distillation theorem which states that with Gaussian operations, it is not 
possible to increase the amount of entanglement (3(| between two parties. Particularly, this 
implies that in each step of an iterative protocol, the bonds need to have at least as much 
entanglement as can be obtained at the output of this step, maximized over all inputs where 
the entanglement is increased. This is indeed a severe restriction, although it does not imply 
the impossibility of such a protocol. One could, e.g., create a hightly entangled state in the 
first step and then approach the desired state by decreasing the entanglement in each step. 
Still, it seems most likely that a sequence of MPS which approach a given state efficiently 
will have to involve more and more bonds simultaneously and thus cannot be constructed 
in an iterative manner. 



D. Correlation functions of Gaussian MPS 

In this section, we show how to compute correlations functions from the maps rM which 
describe the GMPS. We show that this can be done efficiently, i.e., in polynomial time 
independent of the dimension of the graph which is different to the finite dimensional case. 
Of course, this is not too surprising as Gaussian states can be fully characterized by a 
number of paramaters quadratic in the number of modes. 




FIG. 6: If the local operations are described by states r' 1 ' via the Jamiolkowski isomorphism, the 
construction of GMPS can be simplified by replacing the measurement-bond-measurement triples 
by a simple projection onto the EPR state. 
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Let us start with the general case of different as in Fig.[|Ji. The calculation can be 
facilitated by the simple observation that the triples consisting of two projective measure- 
ments and one EPR pair can be replaced by a single projection onto the EPR state, Fig.^b- 
It follows that we can apply the formalism for projective measurements onto the EPR state 
which we presented in Sec. IVIlAl We start from (J^r^. First we partially transpose all 
A modes, then we collapse Ai+i and Bi for all i, and finally we take the Schur complement 
of the merged mode. In case of periodic boundary conditions, this can be expressed by the 
transformation matrix 



n 



1 A K9 B 
t c 



(48) 




which maps ABC onto A'C, where 9b = 0<8 1 is the partial transposition on system B, and 
1Z is the circulant right shift operator, (TZ)%j = Si,j+i mod at <8> 1. Then, the output state, 
i.e., the GMPS characterized by T®, is 



7 = SC A , 



where SCx(U) is the Schur complement of the X part of U, SCx(U) = Uyy — UyxUxxUxy- 
For fixed boundary conditions, the matrix II has to be modified accordingly at the bound- 
aries. All the involved operations scale polynomially in the product NM of the number of 
sites N and the number of modes M. 

In case all the local maps are chosen equal, = T Vi, the above formula can be simplified 
considerably. Therefore, note that the Fourier transform can be taken into the Schur com- 
plement, and that II as well as ©,i =1 rW = r ® l^v are blockwise circulant so that both are 
diagonalized by the Fourier transform. We again adapt the notation of writing the diagonal 
of the Fourier transformed matrices as functions of an angle </>, cf. Sec. 11111 In that case, 
r® 1 is mapped onto the constant function T, and the same holds for 1 and 9 in (|48H . The 
right shift operator 1Z, on the other hand, is transformed to e^t: the EPR measurement 
performed between adjacent sites leads to a complex phase of (j>. Altogether, we have 



n = 



t A e^9 B 
t c 



7 = SC, 



nrn 



Directly expressed in terms of the map T, this reads 

1 



y{4>) -r c -r C | AB A 



ATabiab At 



AB\C 



(49) 



where A = (1^ ; e^Os) is the upper left subblock of ft. 



E. States with rational trigonometric functions as Fourier transforms 

If one restricts to pure MPS (i.e., those for which V is pure) with one mode per site, 
then it follows from Theorem |3 that these states have reflection symmetry, and therefore 
7(0) = 70 + 2 ^2 n>0 7n cos(n0) is real. This implies that the sines in (|49|l can only appear 
in even powers sin 2 ™ = (1 — cos 2 0)". Therefore, the Fourier transform 7 of any pure 
Gaussian MPS, which is a 2 x 2 matrix valued function of 0, has elements which are rational 
functons of cos(^), ("j(4>)) X y = Pxy{cos((f>)) / q X y(cos((f))) with p, q polynomials. The degree of 
the polynomials is limited by the size of Ar^A''', and thus by the number M of the bonds. 
One can easily check that dimp < 2M + 1 and dim q < 2M. 

For the following discussion, let us write those rational functions with a common denom- 
inator d, 

MS\ = 1 ( «( cos ^)) K cos 0)) \ (5Q) 

nV> d(cos(0)) I r(cos(0)) p(cos(0)) I ' 
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where q, p, r, and d are polynomials of degree L. Then, the set of all such 7 with L > 2M+ 1 
encompasses the set of translational invariant GMPS with M bonds. Computing correlation 
functions in a lattice of size N can be done straightforwardly in this representation by taking 
the discrete Fourier transform of "j((f>) which scales polynomially with N, and in the following 
section we show that for one dimension, the correlations can be even computed exactly in 
the limit of an infinite chain. 

It is interesting to note that j(4>) is already determined up to a finite number of possibilities 
by fixing r and d. Since 7 is pure, 1 = det 7 = det 7, and therefore, pq = d 2 +r 2 . Therefore, 
the zeros of pq are the zeros of d 2 + r 2 , such that the only freedom is to choose how to 
distribute the zeros on p and q. On the contrary, fixing only q and d does not give sufficient 
information, while choosing p, q and d (i.e., the diagonal of 7) does not ensure that there 
exists a polynomial r such that pq — r 2 = d 2 . 

From the above, it follows that 2L + 1 parameters are sufficient to describe j(4>), where L 
is still the degree of the polynomials. This encloses all translational invariant Gaussian MPS 
with bond number M < (L — l)/2, which need (2M + l)(2M + 2) = L(L + 1) parameters. 
Therefore, the class of states where 7(0) is a rational function of cos(</>) is a more efficient 
description of translationally invariant states than Gaussian MPS are. 

Let us stress once more that the results of this section hold for arbitrary spatial dimension. 



F. Correlation length 

In the following, we show that the correlations of one-dimensional GMPS decay expo- 
nentially and explicitly derive the correlation length. The derivation only makes use of the 
representation l|50(l of Gaussian MPS and thus holds for the whole class of states where the 
Fourier transform is a rational function of the cosine. We will restrict to the case where the 
state r associated to the GMPS map has only finite entries, which corresponds to the case 
where the denominator d(cos(cf>)) in (|5U[) has no zero on the unit circle. 6 

The correlations are directly obtained by back-transforming the elements of j(4>), which 
are rational functions [7(^)] s = s(cos(0))/d(cos(0)), s € {p, q, r}; in the limit of an infinite 
chain, 



2tt J d(cos(<p)) 



Now transform s, d to complex polynomials via cos(</>) — > (z + 1/ z)/2, and expand with z 



K 



s(z) := z K s(z), d(z) 
in z. Then, 



.,K 



d(z), where K is chosen large enough to make s, d polynomials 



(7s)n = h 



s(z)z r ' 



d(z) 
1 



-dz 

d" 



^ (ui - 1)! dz"*- 1 

Zi-.d(z z )=0 v 



s(z)z 



n-1 



di{z) 



Di 



by the calculus of residues, where z/, is the order of the zero Zi in d and di{z){z— zi) Vi = d{z). 

For n > Ui, Di oc z. v, \ and it follows that the correlations decay exponentially, where the 
correlation length is given by the largest zero of q(z) inside the unit circle. 

This proof holds only for one-dimensional GMPS. However, it can be proven for arbitrary 
spatial dimensions that the correlations decay as o(||ri|| -00 ) by iterated integration by parts 
as in Lemma |5| 



The case where d has zeros on the unit circle corresponds to critical systems, which is why the correlations 
diverge. In the case of a Hamiltonian H = V ffi 1, however, the ground state correlations of P do not 
diverge. As in that case one has pq = d 2 , p/d = d/q need not have a singularity just because q/d has one. 
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G. GMPS as ground states of local Hamiltonians 

Finally, we prove that every Gaussian MPS is the ground state of a local Hamiltonian, 
and show that only a special class of local Hamiltonians has a GMPS as an exact ground 
state. Again, the proof only relies on the representation ffHTTf) - 

Given a state 7 with Fourier transform (|50|) , define H to be the Hamiltonian matrix with 
Fourier transform 

HU) = ( P^ ^)) -r(cos(<f>)) I _ , 51 s 
y -r(cos(0)) g(cos(0)) J 

Then, H corresponds to a local Hamiltonian — the interaction range is the degree of p, q, 1 — 
and according to ©, £{4>) — [y/pq — r 2 ] (cos</>) = (i(cos0), which together with Eq. ijTU|) 
proves that 7 is the ground state of H. 

It is interesting to have a brief look at the converse as well. Given a local Hamiltonian, 
when will it have a GMPS as its ground state? Any local Hamiltonian has a Fourier trans- 
form which consists of polynomials in cos(0), and thus we adapt the notation of Eq. 1)510. 
Then, the ground state is represented by a rational function of cos(^) in Fourier space ex- 
actly if pq — r 2 — d 2 is the square of another polynomial, as can be seen from Eq. I|10|) . In 
terms of the original Hamiltonian, this implies that HqHp — Hq P has to be the square of 
another banded matrix. For example, for the usual case H = V (B 1 one would need V = X 2 
with X again a banded matrix. The Klein-Gordon Hamiltonian (Q, e.g., does not have a 
GMPS as its ground state. 
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APPENDIX A: HAMILTONIAN SIMULATION 



In this Appendix, we discuss the following question. Consider the set S of symplectic 
transformations of N harmonic oscillators on a ring which have both translation and re- 
flection symmetry. Assume that we can implement every local transformation of the form 
Sffi ... 9 S <E S and in addition one element of S corresponding to a nearest-neighbor 
interaction Hamiltonian. Is this set universal for simulating any operation in 5? 



1. Preliminaries 



Lie- Trotter formulae: It is most convenient to discuss the problem on the level of the 
Lie algebra of the symplectic group, since starting with a fixed set of operations in S the 
"reachable" transformations are characterized by the closure of the Lie algebra. That is, if 
e XA and e XB are reachable for every real A then so is e aA+ ^ B and e 7 ^'" 8 ' for all a, (3, 7 € M. 
This follows from the Lie- Trotter formulae 

e *A+pB = Um ( e *A/n e 0B/n\ n and (M) 
n — »oo V / 

e^ = lim ( e A/Vfi e B/^ e -A/Vfi e -B/Vfi) n . (A2 ) 

n — >oo V / 

The Lie algebra: Let S = e tA be a symplectic transformation close to the identity, i.e., 
|<| <C 1. By examining the first order in t we see that the condition SaS T = a is equivalent 
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to Aa = (Aa) T . Hence, we get all possible generators A by multiplying any real symmetric 
matrix with the symplectic matrix a so that 



A=\ „r " W=\ „ %r , (A3) 




where Q, P are symmetric and the block structure corresponds to a direct sum of momentum 
and position space. 

If S G S then Q, P and C have to be circulant matrices (translation symmetry) and 
C = C T (reflection symmetry). Under these symmetries (and utilizing the fact that all 
circulant matrices mutually commute) the blocks of the commutator A" = [A, A'] are given 
by 

P" = 2(C'P-CP'), 

Q" = 2(CQ'-C'Q), (A4) 
C" = Q'P-QP'. 

Note that, even if we do not assume reflection symmetry for A and A' , then A" will have 
it, i.e., C" will be symmetric. For this reason it is crucial to impose reflection symmetry - 
otherwise the set of operations will not be universal. 

For local transformations of the form S® ■ ■ -®S G S the blocks P, Q and C are proportional 
to the identity, since (J),. e Ai — e® > Ai . 

Quadratic Hamiltonians: The relation between the Hamiltonian matrix of a quadratic 
Hamiltonian and the generator of the respective symplectic transformation can be obtained 
from the equation 



I 

Examining the infinitesimal regime yields 

A = aH=( H Q p Hp ). (A6) 

Note that in general a symplectic transformation generated by an A of the form in Eq. (|A3|) 
does not correspond to a semi-bounded Hamiltonian. However, by Eq. (|A1|) we can always 
decompose any symplectic transformation into time evolutions each of which is governed 
by a semi-definite Hamiltonian matrix. Moreover, if we can simulate the time evolution 
governed by a Hamiltonian Tt we can in principle also simulate the evolution according to 
— TL by going to the revival time T rcv for which 

e m{T^-t) = e ~mt_ (A7) 



2. Universality 

Theorem 16 Consider any transformation corresponding to a nearest-neighbor interaction 
Hamiltonian in the set S of translationally invariant symplectic operations with reflection 
symmetry. Together with all local transformations of the form S(B- ■ -®S G S this is universal 
for simulating any transformation in S. 

Proof The theorem is proven in two steps. First we show that using additional local 
transformations only, we can extract each block in Eq. I|A3(I in the sense that one can 
simulate a new operation, where two of the three matrices Q, P, C are set to zero, while the 
third remains unchanged. In this way we can in the second step discuss the reachability for 
each of these three components separately: 

Extracting the blocks. — We will only show how the P-block can be extracted. Similar 
reasoning will then apply for the Q _ block and by utilizing Eq. (|A1(I we can extract the 
C-components by simply subtracting the Q and P blocks. 
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Let us start with a general A of the form in Eq. (|A3|I and first get rid of the C-component. 
To this end consider the commutation relation l|A4(l with C = \\ and Q' = P' = 0. Then 
P" = P, Q" = Q and C" = 0. Starting with an A of this form (C = 0) we may take the 
commutator with Q' = ±1, C = P' = 0. Then P" = Q" = and C" = P. The commutator 
of such a matrix in turn with P' = — il, Q' = C = leads then finally to P" = P with the 
other components zero. 

Generating a basis. — Consider a matrix A which corresponds to a nearest-neighbor in- 
teraction, i.e. in one of the blocks Q,P,C the first off-diagonal is non-zero. By taking 
commutators like in (i) we can easily construct a second generator A' , also corresponding to 
such an interaction, however with the respective off-diagonal in one of the other two blocks. 
By Eq. (|A4I) taking the commutator [A, A'] will lead to a product of the two non-diagonal 
circulant matrices in one of the blocks. This product will again be a circulant matrix, but 
now with a non-zero second off-diagonal. By iterating this procedure we can subsequently 
generate non-zero entries in every off-diagonal and taking linear combinations of these ma- 
trices will thus lead to a basis of symmetric circulant matrices in each of the blocks. Hence, 
together with (i) every generator corresponding to an element of S can be constructed. □ 



3. Remarks 

Revival time and efficiency: Making use of the revival time <| AT|> might be a severe 
handicap concerning the efficiency of the simulation. In particular T rcv of the interaction 
Hamiltonian H mt might grow exponentially with the number N of modes. 

In this case there are two ways to speed up the simulation: either one supplements the set 
of transformations with an additional interaction which is such that it provides an efficient 
simulation of —Hi nt , or one starts with a "good" Hi nt for which this is possible from the 
very beginning. Examples for the latter are the Hamiltonians 

N 

n = +Pf + a(Q l P l+1 +P l Q l+1 ) 1 (A8) 

i=l 
N 

n = 5^(Qi + Q i+ i) 2 + {Pi - P+i) 2 - (A9) 

i=l 

In both cases we can efficiently simulate the evolution according to — H by first applying 
the symplectic transformation Q i— ► P, P i— > —Q and then changing the sign of the diagonal 
in H by local operations. 
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